/* ArithSqrt
 Standard arithmetical integer square root algorithm - G.N.Ward
*/

#include bob:string.format

main()
 {
  local i;
  print( "Standard arithmetical integer square root algorithm\n\n" );
  print( "                                  -----  ArmBob  -----\n" );
  print( "integer N         rounded N      sqrt(N)  rdd(sqrt(N))\n" );
  for (i = 1; i < 31; i ++)
   {
    N = (1 << i) - 1;
    r = arithsqrt( N );
    print( "\nN =",format( N,11,3 ));
    print( "    N =",format( r,6,3 ));
    s = sqrt( 1.0*N ); /* built-in sqrt() */
    t = floor( s + 0.5 );
    print(format( s,12,2 ));
    print(format( t,8,3 ));
   }
   print( "\n" );
 }
 
arithsqrt(b)
 {
   u = b;
   v = 1 << 6;
   n = 3;
   repeat
    {
      u >>= 6;
      if (u != 0)
       { n += 3; v <<= 6; }
    }
   until (u == 0);
   if (n > 15)
    { v = 1 << 30; n = 15; }
   v = 1 << 30; n = 15;
   repeat
    {
      d = u + v;
      u >>= 1;
      if (b >= d) { b -= d; u += v; }
      v >>= 2;
      n --;
    }
   until (n < 0);
   if (b > u) u ++;
   return u;
 }
     
format(n,l,p)
{
  local dec,i,int,neg,s;
  s = "";
  switch (typeof(n))
    {
      case REAL:
        if (neg = (n < 0.0)) n = -n;
        int = floor(n);
        if (neg) s = "-";
        s += int_10(int);
        s += ".";
        dec = n - 1.0*int; 
        for (i = 0; i < p; i ++) dec *= 10.0;
        dec = floor(dec + 0.5);
        d = "";
        for (i = 0;i < p; i ++)
          { d  = '0' + dec%10 + d; dec /= 10; }
        s += d;
        break;
      case INTEGER:
        if (neg = (n < 0)) n = -n;
        if (neg) s = "-";
        s += int_10(n);
        break;
      case STRING:
        s = n;        
        break;
      default:
        s = "first argument not real, integer or string";        
    }
  for (i = l - sizeof(s); i > 0; i --) s = " " + s;
  return s;
}
 
int_10(n)
{
  local s;
  s = "";
  if (n == 0) s = "0";
  else while(n) { s = '0' + n%10 + s; n /= 10; }
  return s;
}

